[Last 2000909, Final 191205-1212, 18project, 3T3] MakeCountmatrix_forGEO_201012 TSSpm5kb__3T3_EGFP18_Dox_20200811__mainH3p3_ChIL01100111_ATAC0049

Setup

library(ggplot2)
library(dplyr)
library(tidyr)
library(purrr)
source("/home/guestA/n70275b/work/rscripts/geomNorm.R")

# Helper function
#ggpoints <- function(x,...) 
#  ggplot(x,...) + geom_point(size=3,stroke=1) +
#  ggrepel::geom_text_repel(size=4) + theme_minimal() + mycolor

## ラベルあり
ggpoints <- function(x,...) 
  ggplot(x,...) + geom_point(stroke=1) +
  ggrepel::geom_text_repel(size=4) + theme_minimal() + mycolor

## ラベルなし
#ggpoints <- function(x,...) 
#  ggplot(x,...) + geom_point(stroke=1) + theme_minimal() + mycolor


print(Sys.Date())

[1] "2020-10-13"

print(sessionInfo(),locale=FALSE)

R version 4.0.1 (2020-06-06)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux

Matrix products: default
BLAS/LAPACK: /usr/local/intel2018_up1/compilers_and_libraries_2018.0.128/linux/mkl/lib/intel64_lin/libmkl_intel_lp64.so

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] purrr_0.3.4   tidyr_1.1.2   ggplot2_3.3.2 dplyr_1.0.2  

loaded via a namespace (and not attached):
  [1] colorspace_1.4-1            ellipsis_0.3.1              XVector_0.28.0              GenomicRanges_1.40.0        rstudioapi_0.11             farver_2.0.3                ggrepel_0.8.2              
  [8] DT_0.15                     bit64_4.0.5                 AnnotationDbi_1.50.3        fansi_0.4.1                 xml2_1.3.2                  R.methodsS3_1.8.1           splines_4.0.1              
 [15] geneplotter_1.66.0          knitr_1.30                  jsonlite_1.7.1              Rsamtools_2.4.0             seqLogo_1.54.3              annotate_1.66.0             GO.db_3.11.4               
 [22] dbplyr_1.4.4                png_0.1-7                   R.oo_1.24.0                 pheatmap_1.0.12             uwot_0.1.8                  shiny_1.5.0                 msigdbr_7.2.1              
 [29] readr_1.4.0                 compiler_4.0.1              httr_1.4.2                  lazyeval_0.2.2              assertthat_0.2.1            Matrix_1.2-18               fastmap_1.0.1              
 [36] cli_2.0.2                   later_1.1.0.1               htmltools_0.5.0             prettyunits_1.1.1           tools_4.0.1                 TFMPvalue_0.0.8             gtable_0.3.0               
 [43] glue_1.4.2                  GenomeInfoDbData_1.2.3      reshape2_1.4.4              rappdirs_0.3.1              Rcpp_1.0.5                  Biobase_2.48.0              vctrs_0.3.4                
 [50] Biostrings_2.56.0           rtracklayer_1.48.0          CNEr_1.24.0                 xfun_0.18                   stringr_1.4.0               mime_0.9                    miniUI_0.1.1.1             
 [57] lifecycle_0.2.0             poweRlaw_0.70.6             gtools_3.8.2                XML_3.99-0.5                zlibbioc_1.34.0             scales_1.1.1                BSgenome_1.56.0            
 [64] hms_0.5.3                   promises_1.1.1              parallel_4.0.1              SummarizedExperiment_1.18.2 RColorBrewer_1.1-2          yaml_2.2.1                  curl_4.3                   
 [71] memoise_1.1.0               gridExtra_2.3               chromVAR_1.10.0             biomaRt_2.45.5              stringi_1.5.3               RSQLite_2.2.1               genefilter_1.70.0          
 [78] S4Vectors_0.26.1            caTools_1.18.0              BiocGenerics_0.34.0         BiocParallel_1.22.0         GenomeInfoDb_1.24.2         rlang_0.4.8                 pkgconfig_2.0.3            
 [85] matrixStats_0.57.0          bitops_1.0-6                pracma_2.2.9                evaluate_0.14               lattice_0.20-41             GenomicAlignments_1.24.0    htmlwidgets_1.5.2          
 [92] labeling_0.3                bit_4.0.4                   tidyselect_1.1.0            ggsci_2.9                   plyr_1.8.6                  magrittr_1.5                DESeq2_1.28.1              
 [99] R6_2.4.1                    IRanges_2.22.2              generics_0.0.2              DelayedArray_0.14.1         DBI_1.1.0                   pillar_1.4.6                withr_2.3.0                
[106] KEGGREST_1.28.0             survival_3.2-7              RCurl_1.98-1.2              tibble_3.0.3                crayon_1.3.4                utf8_1.1.4                  plotly_4.9.2.1             
[113] BiocFileCache_1.12.1        rmarkdown_2.4               viridis_0.5.1               progress_1.2.2              TFBSTools_1.26.0            locfit_1.5-9.4              grid_4.0.1                 
[120] data.table_1.13.0           blob_1.2.1                  digest_0.6.25               xtable_1.8-4                httpuv_1.5.4                R.utils_2.10.1              DirichletMultinomial_1.30.0
[127] openssl_1.4.3               stats4_4.0.1                munsell_0.5.0               viridisLite_0.3.0           askpass_1.1                

select <- dplyr::select
rename <- dplyr::rename #191203
count <- dplyr::count #191203

Parameters

modify here

# Files


#deftable <- "~/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Rserver_191203/deftable_BRB_noumi_new_190520_fin191205ver.txt" #Umi補正なし (BRB)

deftable <- "/home/guestA/o70578a/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Last_Rserver_200811/deftable_BRB_noumi_new_190520_Last20200811ver.txt"

#deftable <- "~/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Rserver_191203/deftable_BRB_noumi_new_190520_fin191205ver.txt" #Umi補正なし (BRB)


#deftable <- "deftable_BRB_noumi_new_190520.txt" #Umi補正なし (BRB)

#deftable <- "~/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/deftable_BRB_noumi_new_190520.txt" #Umi補正なし (BRB)


## Data selection (filter rows of deftable)
#use <- quo(!grepl("^18",group) & (group != "Nc-minusTryd"))
#use <- quo(TRUE) # use all
use <- quo(type != "C2C12")

# Species specific parameters
species <- "Mus musculus"
biomartann <- "mmusculus_gene_ensembl"
maxchrom <- 19 # 19: mouse, 22: human


# Graphics
# aesthetic mapping of labels
#myaes <- aes(colour=enzyme,shape=leg,label=rep) 
#myaes <- aes(colour=growth,shape=type,size=count) #ラベルなし
#myaes <- aes(colour=growth,shape=type,label=replicate,size=count) #ラベルあり
#myaes <- aes(colour=enzyme,shape=leg,label=replicate) #ラベルあり
#myaes <- aes(colour=enzyme,shape=leg,label=factor(rep))
#myaes <- aes(colour=type, shape=revcro, label=read, size=count)
#myaes <- aes(colour=type, shape=revcro, label=read)

#myaes <- aes(colour=growth,shape=type,label=replicate,size=count) #ラベルあり
#myaes <- aes(colour=time,shape=type,size=count,label=replicate)
#myaes <- aes(colour=WT_KO_intact_CTX, shape=Day,size=count,label=f_m)

#myaes <- aes(colour=WT_KO_intact_CTX, shape=Day, label=f_m) #サイズを変えず
#myaes <- aes(colour=growth,shape=type,label=replicate,size=count) #ラベルあり
myaes <- aes(colour=time,shape=type,label=rep,size=count) #ラベルあり
myaes2 <- aes(colour=time,shape=type) #kuwa add

# color palette of points: See vignette("ggsci")
mycolor <- ggsci::scale_color_aaas()



# PCA/UMAP
scalerows <- TRUE # gene-wise scaling (pattern is the matter?)
ntop <- 500 # number of top-n genes with high variance
seed <- 123 # set another number if UMAP looks not good
n_nei <- 6  # number of neighboring data points in UMAP #ここをどうしたらいい？


# DESeq2
#model <- ~groupn+lead #dateも追加
#model <- ~leg + enzyme + leg:enzyme
#model <- ~type+growth#+type:growth
#model <- ~group+lead


#model <- ~group
#model <- ~type+growth+type:growth #これでは相互作用が入っていない
#model <- ~type+growth #これでは相互作用が入っていない


model <- ~group
#model <- ~type+growth+growth:type

fdr <- 0.1 # acceptable false discovery rate
lfcthreth <- log2(1) # threshold in abs(log2FC)

# controls should be placed in the right
contrast <- list(
  
  group_UI_Doxplus_vs_minus = c("group", "BRB_UI_DoxPlus", "BRB_UI_DoxMinus"),
  group_0h_Doxplus_vs_minus = c("group", "BRB_0h_DoxPlus", "BRB_0h_DoxMinus"),
  group_24h_Doxplus_vs_minus = c("group", "BRB_24h_DoxPlus", "BRB_24h_DoxMinus"),
  group_48h_Doxplus_vs_minus = c("group", "BRB_48h_DoxPlus", "BRB_48h_DoxMinus")
  
  
  #group_UI_Doxplus_vs_minus = c("group", "Doxplus_UI", "Doxminus_UI"),
  #group_Diff0h_Doxplus_vs_minus = c("group", "Doxplus_Diff0h", "Doxminus_Diff0h"),
  #group_Diff24h_Doxplus_vs_minus = c("group", "Doxplus_Diff24h", "Doxminus_Diff24h"),
  #group_Diff48h_Doxplus_vs_minus = c("group", "Doxplus_Diff48h", "Doxminus_Diff48h")
  
  
  #Intercept = list("Intercept"), # reference level
  #leg_LvsR = c("leg", "L", "R"),
  #enz_KvsC = c("enzyme","K","C")
  #legL.enzK = list("legL.enzymeK") # interaction
  
  #type_Doxplus_vs_minus = c("type", "Doxplus", "Doxminus")
)

Retrieve Biomart

#if(!exists("e2g")){
#  #ensembl <- biomaRt::useMart("ENSEMBL_MART_ENSEMBL",host="asia.ensembl.org")
#  #ensembl <- biomaRt::useMart("ENSEMBL_MART_ENSEMBL",host="uswest.ensembl.org")
#  ensembl <- biomaRt::useMart("ENSEMBL_MART_ENSEMBL",host="useast.ensembl.org")
#  mart <- biomaRt::useDataset(biomartann,mart=ensembl)
#  e2g <- biomaRt::getBM(attributes=c("ensembl_gene_id","external_gene_name",
#    "gene_biotype","chromosome_name"), mart=mart) %>% as_tibble %>%
#  rename(
#    ens_gene = ensembl_gene_id,
#    ext_gene = external_gene_name,
#    biotype = gene_biotype,
#    chr = chromosome_name
#  )
#}

e2g <- readr::read_csv("/home/guestA/o70578a/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Last_Rserver_200811/ensemble_list_useast.csv")

─ Column specification ──────────────────────────────────────────────────────────────────────────────────────────────
cols(
  ens_gene = col_character(),
  ext_gene = col_character(),
  biotype = col_character(),
  chr = col_character()
)

annotate <- partial(right_join,e2g,by="ens_gene")

#annotate <- partial(right_join,e2g,by="ens_gene")

#-----#
nrow(e2g)

[1] 56305

#readr::write_csv(e2g,"ensemble_list_asia.csv")
#readr::write_csv(e2g,"ensemble_list_uswest.csv")
#readr::write_csv(e2g,"ensemble_list_useast.csv")

Load counts

filepath_summary <- "/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/ChIL01100111_ATAC0049L1__3T3_EGFP18_Dox__TSS_pm5kb_20200624.count.txt"

filepath_selectedTSScount <- "/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/ChIL01100111_ATAC0049L1__3T3_EGFP18_Dox__TSS_pm5kb_20200624.count__count_firstTSS_select.csv"


colnames_title <- c("ens_gene","chr","TSSstart","TSSend","strand","TSS","Start","End")
  
samples_ATAC <- c(
  "ATAC_UI_DoxMinus_1","ATAC_UI_DoxMinus_2","ATAC_UI_DoxMinus_3","ATAC_UI_DoxMinus_4",
  "ATAC_UI_DoxPlus_1","ATAC_UI_DoxPlus_2","ATAC_UI_DoxPlus_3","ATAC_UI_DoxPlus_4",
  "ATAC_48h_DoxMinus_1","ATAC_48h_DoxMinus_2","ATAC_48h_DoxMinus_3","ATAC_48h_DoxMinus_4",
  "ATAC_48h_DoxPlus_1","ATAC_48h_DoxPlus_2","ATAC_48h_DoxPlus_3","ATAC_48h_DoxPlus_4")

samples_H3p3 <- c(
  "H3p3_UI_DoxMinus_1","H3p3_UI_DoxMinus_2","H3p3_UI_DoxPlus_1","H3p3_UI_DoxPlus_2",
  "H3p3_0h_DoxMinus_1","H3p3_0h_DoxMinus_2","H3p3_0h_DoxPlus_1","H3p3_0h_DoxPlus_2",
  "H3p3_24h_DoxMinus_1","H3p3_24h_DoxMinus_2","H3p3_24h_DoxPlus_1","H3p3_24h_DoxPlus_2",
  "H3p3_48h_DoxMinus_1","H3p3_48h_DoxMinus_2","H3p3_48h_DoxPlus_1","H3p3_48h_DoxPlus_2")

samples_H3K27ac <- c(
  "H3K27ac_UI_DoxMinus_1","H3K27ac_UI_DoxMinus_2","H3K27ac_UI_DoxPlus_1","H3K27ac_UI_DoxPlus_2",
  "H3K27ac_0h_DoxMinus_1","H3K27ac_0h_DoxMinus_2","H3K27ac_0h_DoxPlus_1","H3K27ac_0h_DoxPlus_2",
  "H3K27ac_24h_DoxMinus_1","H3K27ac_24h_DoxMinus_2","H3K27ac_24h_DoxPlus_1","H3K27ac_24h_DoxPlus_2",
  "H3K27ac_48h_DoxMinus_1","H3K27ac_48h_DoxMinus_2","H3K27ac_48h_DoxPlus_1","H3K27ac_48h_DoxPlus_2")

samples_H3K4me3 <- c(
  "H3K4me3_UI_DoxMinus_1","H3K4me3_UI_DoxMinus_2","H3K4me3_UI_DoxPlus_1","H3K4me3_UI_DoxPlus_2",
  "H3K4me3_0h_DoxMinus_1","H3K4me3_0h_DoxMinus_2","H3K4me3_0h_DoxPlus_1","H3K4me3_0h_DoxPlus_2",
  "H3K4me3_24h_DoxMinus_1","H3K4me3_24h_DoxMinus_2","H3K4me3_24h_DoxPlus_1","H3K4me3_24h_DoxPlus_2",
  "H3K4me3_48h_DoxMinus_1","H3K4me3_48h_DoxMinus_2","H3K4me3_48h_DoxPlus_1","H3K4me3_48h_DoxPlus_2")


selectedTSScount <- readr::read_csv(filepath_selectedTSScount) %>% mutate_if(is.double, as.integer)

─ Column specification ──────────────────────────────────────────────────────────────────────────────────────────────
cols(
  .default = col_double(),
  chr = col_character(),
  ens_gene = col_character(),
  score = col_character(),
  strand = col_character()
)
ℹ Use `spec()` for the full column specifications.

selectedTSScount1 <- selectedTSScount %>% dplyr::select(all_of(colnames_title), all_of(samples_H3p3), all_of(samples_H3K4me3),all_of(samples_H3K27ac), all_of(samples_ATAC)) %>% dplyr::rename(TSSminus5kb=TSSstart,TSSplus5kb=TSSend)

nrow(selectedTSScount)

[1] 55450

colnames(selectedTSScount)

 [1] "chr"                     "ens_gene"                "TSSstart"                "TSSend"                  "score"                   "strand"                  "TSS"                    
 [8] "Start"                   "End"                     "position"                "ATAC_UI_DoxMinus_1"      "ATAC_UI_DoxMinus_2"      "ATAC_UI_DoxMinus_3"      "ATAC_UI_DoxMinus_4"     
[15] "ATAC_UI_DoxPlus_1"       "ATAC_UI_DoxPlus_2"       "ATAC_UI_DoxPlus_3"       "ATAC_UI_DoxPlus_4"       "ATAC_48h_DoxMinus_1"     "ATAC_48h_DoxMinus_2"     "ATAC_48h_DoxMinus_3"    
[22] "ATAC_48h_DoxMinus_4"     "ATAC_48h_DoxPlus_1"      "ATAC_48h_DoxPlus_2"      "ATAC_48h_DoxPlus_3"      "ATAC_48h_DoxPlus_4"      "H3p3_UI_DoxMinus_1"      "H3p3_UI_DoxMinus_2"     
[29] "H3p3_UI_DoxPlus_1"       "H3p3_UI_DoxPlus_2"       "H3p3_0h_DoxMinus_1"      "H3p3_0h_DoxMinus_2"      "H3p3_0h_DoxPlus_1"       "H3p3_0h_DoxPlus_2"       "H3p3_24h_DoxMinus_1"    
[36] "H3p3_24h_DoxMinus_2"     "H3p3_24h_DoxPlus_1"      "H3p3_24h_DoxPlus_2"      "H3p3_48h_DoxMinus_1"     "H3p3_48h_DoxMinus_2"     "H3p3_48h_DoxPlus_1"      "H3p3_48h_DoxPlus_2"     
[43] "H3K27ac_UI_DoxMinus_1"   "H3K27ac_UI_DoxMinus_2"   "H3K27ac_UI_DoxPlus_1"    "H3K27ac_UI_DoxPlus_2"    "H3K27ac_0h_DoxMinus_1"   "H3K27ac_0h_DoxMinus_2"   "H3K27ac_0h_DoxPlus_1"   
[50] "H3K27ac_0h_DoxPlus_2"    "H3K27ac_24h_DoxMinus_1"  "H3K27ac_24h_DoxMinus_2"  "H3K27ac_24h_DoxPlus_1"   "H3K27ac_24h_DoxPlus_2"   "H3K27ac_48h_DoxMinus_1"  "H3K27ac_48h_DoxMinus_2" 
[57] "H3K27ac_48h_DoxPlus_1"   "H3K27ac_48h_DoxPlus_2"   "H3K4me3_UI_DoxMinus_1"   "H3K4me3_UI_DoxMinus_2"   "H3K4me3_UI_DoxPlus_1"    "H3K4me3_UI_DoxPlus_2"    "H3K4me3_0h_DoxMinus_1"  
[64] "H3K4me3_0h_DoxMinus_2"   "H3K4me3_0h_DoxPlus_1"    "H3K4me3_0h_DoxPlus_2"    "H3K4me3_24h_DoxMinus_1"  "H3K4me3_24h_DoxMinus_2"  "H3K4me3_24h_DoxPlus_1"   "H3K4me3_24h_DoxPlus_2"  
[71] "H3K4me3_48h_DoxMinus_1"  "H3K4me3_48h_DoxMinus_2"  "H3K4me3_48h_DoxPlus_1"   "H3K4me3_48h_DoxPlus_2"   "H3K27me3_UI_DoxMinus_1"  "H3K27me3_UI_DoxMinus_2"  "H3K27me3_UI_DoxPlus_1"  
[78] "H3K27me3_UI_DoxPlus_2"   "H3K27me3_0h_DoxMinus_1"  "H3K27me3_0h_DoxMinus_2"  "H3K27me3_0h_DoxPlus_1"   "H3K27me3_0h_DoxPlus_2"   "H3K27me3_24h_DoxMinus_1" "H3K27me3_24h_DoxMinus_2"
[85] "H3K27me3_24h_DoxPlus_1"  "H3K27me3_24h_DoxPlus_2"  "H3K27me3_48h_DoxMinus_1" "H3K27me3_48h_DoxMinus_2" "H3K27me3_48h_DoxPlus_1"  "H3K27me3_48h_DoxPlus_2" 

####
nrow(selectedTSScount1)

[1] 55450

colnames(selectedTSScount1)

 [1] "ens_gene"               "chr"                    "TSSminus5kb"            "TSSplus5kb"             "strand"                 "TSS"                    "Start"                  "End"                   
 [9] "H3p3_UI_DoxMinus_1"     "H3p3_UI_DoxMinus_2"     "H3p3_UI_DoxPlus_1"      "H3p3_UI_DoxPlus_2"      "H3p3_0h_DoxMinus_1"     "H3p3_0h_DoxMinus_2"     "H3p3_0h_DoxPlus_1"      "H3p3_0h_DoxPlus_2"     
[17] "H3p3_24h_DoxMinus_1"    "H3p3_24h_DoxMinus_2"    "H3p3_24h_DoxPlus_1"     "H3p3_24h_DoxPlus_2"     "H3p3_48h_DoxMinus_1"    "H3p3_48h_DoxMinus_2"    "H3p3_48h_DoxPlus_1"     "H3p3_48h_DoxPlus_2"    
[25] "H3K4me3_UI_DoxMinus_1"  "H3K4me3_UI_DoxMinus_2"  "H3K4me3_UI_DoxPlus_1"   "H3K4me3_UI_DoxPlus_2"   "H3K4me3_0h_DoxMinus_1"  "H3K4me3_0h_DoxMinus_2"  "H3K4me3_0h_DoxPlus_1"   "H3K4me3_0h_DoxPlus_2"  
[33] "H3K4me3_24h_DoxMinus_1" "H3K4me3_24h_DoxMinus_2" "H3K4me3_24h_DoxPlus_1"  "H3K4me3_24h_DoxPlus_2"  "H3K4me3_48h_DoxMinus_1" "H3K4me3_48h_DoxMinus_2" "H3K4me3_48h_DoxPlus_1"  "H3K4me3_48h_DoxPlus_2" 
[41] "H3K27ac_UI_DoxMinus_1"  "H3K27ac_UI_DoxMinus_2"  "H3K27ac_UI_DoxPlus_1"   "H3K27ac_UI_DoxPlus_2"   "H3K27ac_0h_DoxMinus_1"  "H3K27ac_0h_DoxMinus_2"  "H3K27ac_0h_DoxPlus_1"   "H3K27ac_0h_DoxPlus_2"  
[49] "H3K27ac_24h_DoxMinus_1" "H3K27ac_24h_DoxMinus_2" "H3K27ac_24h_DoxPlus_1"  "H3K27ac_24h_DoxPlus_2"  "H3K27ac_48h_DoxMinus_1" "H3K27ac_48h_DoxMinus_2" "H3K27ac_48h_DoxPlus_1"  "H3K27ac_48h_DoxPlus_2" 
[57] "ATAC_UI_DoxMinus_1"     "ATAC_UI_DoxMinus_2"     "ATAC_UI_DoxMinus_3"     "ATAC_UI_DoxMinus_4"     "ATAC_UI_DoxPlus_1"      "ATAC_UI_DoxPlus_2"      "ATAC_UI_DoxPlus_3"      "ATAC_UI_DoxPlus_4"     
[65] "ATAC_48h_DoxMinus_1"    "ATAC_48h_DoxMinus_2"    "ATAC_48h_DoxMinus_3"    "ATAC_48h_DoxMinus_4"    "ATAC_48h_DoxPlus_1"     "ATAC_48h_DoxPlus_2"     "ATAC_48h_DoxPlus_3"     "ATAC_48h_DoxPlus_4"    

print(selectedTSScount1)

#### save
TSScount_writefile1 <- basename(filepath_selectedTSScount) %>% gsub("ChIL01100111","ChILseq",.) %>% gsub("ATAC0049L1","ATACseq",.) %>% gsub(".count__count_firstTSS_select","__CountMatrix",.) %>% gsub(".count__count_firstTSS_select","___CountMatrix",.)
print(TSScount_writefile1)

[1] "ChILseq_ATACseq__3T3_EGFP18_Dox__TSS_pm5kb_20200624__CountMatrix.csv"

selectedTSScount1 %>% readr::write_csv(TSScount_writefile1) ##これを登録
##

#それぞれのlog2FC (読み込みが大変なため、そのまま登録)
filepath_log2FC_H3p3 <- "/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/2gun/TSS_pm5kb_20200624_resultsall_fdr0p1_H3p3.csv"

filepath_log2FC_H3K4me3 <- "/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/2gun/TSS_pm5kb_20200624_resultsall_fdr0p1_H3K4me3.csv"

filepath_log2FC_H3K27ac <- "/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/2gun/TSS_pm5kb_20200624_resultsall_fdr0p1_H3K27ac.csv"

filepath_log2FC_ATAC <- "/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/2gun/TSS_pm5kb_20200624_resultsall_fdr0p1_ATAC.csv"

####

log2FC_H3p3 <- readr::read_csv(filepath_log2FC_H3p3)

─ Column specification ──────────────────────────────────────────────────────────────────────────────────────────────
cols(
  aspect = col_character(),
  ens_gene = col_character(),
  ext_gene = col_character(),
  biotype = col_character(),
  chr = col_character(),
  baseMean = col_double(),
  log2FoldChange = col_double(),
  lfcSE = col_double(),
  stat = col_double(),
  pvalue = col_double(),
  padj = col_double()
)

colnames(log2FC_H3p3)

 [1] "aspect"         "ens_gene"       "ext_gene"       "biotype"        "chr"            "baseMean"       "log2FoldChange" "lfcSE"          "stat"           "pvalue"         "padj"          

log2FC_H3p3 %>% group_by(aspect) %>%  summarise(n())

`summarise()` ungrouping output (override with `.groups` argument)

####
log2FC_H3K4me3 <- readr::read_csv(filepath_log2FC_H3K4me3)

─ Column specification ──────────────────────────────────────────────────────────────────────────────────────────────
cols(
  aspect = col_character(),
  ens_gene = col_character(),
  ext_gene = col_character(),
  biotype = col_character(),
  chr = col_character(),
  baseMean = col_double(),
  log2FoldChange = col_double(),
  lfcSE = col_double(),
  stat = col_double(),
  pvalue = col_double(),
  padj = col_double()
)

colnames(log2FC_H3K4me3)

 [1] "aspect"         "ens_gene"       "ext_gene"       "biotype"        "chr"            "baseMean"       "log2FoldChange" "lfcSE"          "stat"           "pvalue"         "padj"          

log2FC_H3K4me3 %>% group_by(aspect) %>%  summarise(n())

`summarise()` ungrouping output (override with `.groups` argument)

####
log2FC_H3K27ac <- readr::read_csv(filepath_log2FC_H3K27ac)

─ Column specification ──────────────────────────────────────────────────────────────────────────────────────────────
cols(
  aspect = col_character(),
  ens_gene = col_character(),
  ext_gene = col_character(),
  biotype = col_character(),
  chr = col_character(),
  baseMean = col_double(),
  log2FoldChange = col_double(),
  lfcSE = col_double(),
  stat = col_double(),
  pvalue = col_double(),
  padj = col_double()
)

colnames(log2FC_H3K27ac)

 [1] "aspect"         "ens_gene"       "ext_gene"       "biotype"        "chr"            "baseMean"       "log2FoldChange" "lfcSE"          "stat"           "pvalue"         "padj"          

log2FC_H3K27ac %>% group_by(aspect) %>%  summarise(n())

`summarise()` ungrouping output (override with `.groups` argument)

####
log2FC_ATAC <- readr::read_csv(filepath_log2FC_ATAC)

─ Column specification ──────────────────────────────────────────────────────────────────────────────────────────────
cols(
  aspect = col_character(),
  ens_gene = col_character(),
  ext_gene = col_character(),
  biotype = col_character(),
  chr = col_character(),
  baseMean = col_double(),
  log2FoldChange = col_double(),
  lfcSE = col_double(),
  stat = col_double(),
  pvalue = col_double(),
  padj = col_double()
)

colnames(log2FC_ATAC)

 [1] "aspect"         "ens_gene"       "ext_gene"       "biotype"        "chr"            "baseMean"       "log2FoldChange" "lfcSE"          "stat"           "pvalue"         "padj"          

log2FC_ATAC %>% group_by(aspect) %>%  summarise(n())

`summarise()` ungrouping output (override with `.groups` argument)

####

log2FC_correlationlist <- readr::read_csv("/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/log2FC/tables/Spread_log2FC_ChILATACBRB_cutoff10__withCluster__H3p3clus3.csv")

─ Column specification ──────────────────────────────────────────────────────────────────────────────────────────────
cols(
  ens_gene = col_character(),
  ext_gene = col_character(),
  biotype = col_character(),
  chr = col_character(),
  time = col_character(),
  H3p3 = col_double(),
  H3K4me3 = col_double(),
  H3K27ac = col_double(),
  H3K27me3 = col_double(),
  ATAC = col_double(),
  BRB = col_double(),
  BRBgroupMean = col_double(),
  H3p3cluster = col_double(),
  BRBDEGcluster = col_double()
)

log2FC_corrtable <- log2FC_correlationlist %>% dplyr::select(-H3K27me3) %>% dplyr::rename(H3p3_log2FC_Doxplus_vs_minus=H3p3,  H3K4me3_log2FC_Doxplus_vs_minus=H3K4me3, H3K27ac_log2FC_Doxplus_vs_minus=H3K27ac, ATAC_log2FC_Doxplus_vs_minus=ATAC, BRB_log2FC_Doxplus_vs_minus=BRB) %>% dplyr::rename(BRB_normalizedcount_groupMean=BRBgroupMean)

nrow(log2FC_correlationlist)

[1] 13452

colnames(log2FC_correlationlist)

 [1] "ens_gene"      "ext_gene"      "biotype"       "chr"           "time"          "H3p3"          "H3K4me3"       "H3K27ac"       "H3K27me3"      "ATAC"          "BRB"           "BRBgroupMean" 
[13] "H3p3cluster"   "BRBDEGcluster"

nrow(log2FC_corrtable)

[1] 13452

colnames(log2FC_corrtable)

 [1] "ens_gene"                        "ext_gene"                        "biotype"                         "chr"                             "time"                            "H3p3_log2FC_Doxplus_vs_minus"   
 [7] "H3K4me3_log2FC_Doxplus_vs_minus" "H3K27ac_log2FC_Doxplus_vs_minus" "ATAC_log2FC_Doxplus_vs_minus"    "BRB_log2FC_Doxplus_vs_minus"     "BRB_normalizedcount_groupMean"   "H3p3cluster"                    
[13] "BRBDEGcluster"                  

log2FC_corrtable %>% group_by(time) %>% summarise(count=n())

`summarise()` ungrouping output (override with `.groups` argument)

log2FC_corrtable %>% filter(!is.na(BRB_normalizedcount_groupMean)) %>% group_by(time) %>% summarise(count=n())

`summarise()` ungrouping output (override with `.groups` argument)

###
log2FC_corrtable %>% group_by(time) %>% filter(BRBDEGcluster=="3") %>% summarise(count=n())

`summarise()` ungrouping output (override with `.groups` argument)

log2FC_corrtable %>% filter(!is.na(BRB_normalizedcount_groupMean)) %>% group_by(time) %>% filter(BRBDEGcluster=="3") %>% summarise(count=n())

`summarise()` ungrouping output (override with `.groups` argument)

###
log2FC_corrtable %>% group_by(time) %>% filter(BRBDEGcluster=="3") %>% filter(ext_gene=="Acta1")

NA

"/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__zscore_All.csv"

[1] "/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__zscore_All.csv"

#それぞれのNormCount (読み込みが大変なため、そのまま登録)
filepath_Normcount_H3p3 <- "/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_H3p3_genename.csv"
filepath_Normcount_H3K4me3 <- "/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_H3K4me3_genename.csv"
filepath_Normcount_H3K27ac <- "/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_H3K27ac_genename.csv"
filepath_Normcount_ATAC <- "/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_ATAC_genename.csv"

readr::read_csv(filepath_Normcount_H3p3)

─ Column specification ──────────────────────────────────────────────────────────────────────────────────────────────
cols(
  .default = col_double(),
  ens_gene = col_character(),
  ext_gene = col_character(),
  biotype = col_character()
)
ℹ Use `spec()` for the full column specifications.

4182 parsing failures.
  row col expected actual                                                                                                                                                                                                                                        file
51016 chr a double      X '/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_H3p3_genename.csv'
51017 chr a double      X '/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_H3p3_genename.csv'
51018 chr a double      X '/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_H3p3_genename.csv'
51019 chr a double      X '/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_H3p3_genename.csv'
51020 chr a double      X '/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_H3p3_genename.csv'
..... ... ........ ...... ...........................................................................................................................................................................................................................................
See problems(...) for more details.

readr::read_csv(filepath_Normcount_H3K4me3)

─ Column specification ──────────────────────────────────────────────────────────────────────────────────────────────
cols(
  .default = col_double(),
  ens_gene = col_character(),
  ext_gene = col_character(),
  biotype = col_character()
)
ℹ Use `spec()` for the full column specifications.

4182 parsing failures.
  row col expected actual                                                                                                                                                                                                                                           file
51016 chr a double      X '/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_H3K4me3_genename.csv'
51017 chr a double      X '/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_H3K4me3_genename.csv'
51018 chr a double      X '/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_H3K4me3_genename.csv'
51019 chr a double      X '/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_H3K4me3_genename.csv'
51020 chr a double      X '/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_H3K4me3_genename.csv'
..... ... ........ ...... ..............................................................................................................................................................................................................................................
See problems(...) for more details.

readr::read_csv(filepath_Normcount_H3K27ac)

─ Column specification ──────────────────────────────────────────────────────────────────────────────────────────────
cols(
  .default = col_double(),
  ens_gene = col_character(),
  ext_gene = col_character(),
  biotype = col_character()
)
ℹ Use `spec()` for the full column specifications.

4182 parsing failures.
  row col expected actual                                                                                                                                                                                                                                           file
51016 chr a double      X '/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_H3K27ac_genename.csv'
51017 chr a double      X '/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_H3K27ac_genename.csv'
51018 chr a double      X '/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_H3K27ac_genename.csv'
51019 chr a double      X '/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_H3K27ac_genename.csv'
51020 chr a double      X '/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_H3K27ac_genename.csv'
..... ... ........ ...... ..............................................................................................................................................................................................................................................
See problems(...) for more details.

readr::read_csv(filepath_Normcount_ATAC)

─ Column specification ──────────────────────────────────────────────────────────────────────────────────────────────
cols(
  .default = col_double(),
  ens_gene = col_character(),
  ext_gene = col_character(),
  biotype = col_character()
)
ℹ Use `spec()` for the full column specifications.

4182 parsing failures.
  row col expected actual                                                                                                                                                                                                                                        file
51016 chr a double      X '/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_ATAC_genename.csv'
51017 chr a double      X '/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_ATAC_genename.csv'
51018 chr a double      X '/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_ATAC_genename.csv'
51019 chr a double      X '/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_ATAC_genename.csv'
51020 chr a double      X '/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/R_server_Last_20200811_BRB_ATAC_ChIL_fistExon_mainH3p3/TSS_pm5kb_20200624__normcount_ATAC_genename.csv'
..... ... ........ ...... ...........................................................................................................................................................................................................................................
See problems(...) for more details.